clear all;
clc;
tstart=1; %start time
tmax=50; %max time of simulation
time=tstart:tmax; %time vector

n_iterations=1;
%nparticles = 1000;
nparticles = 10;
npos=2*length(time)+3; %size of pos arrays = 2t+3
pos=zeros(npos,npos,npos);
pos_x=1:npos;
pos_y=1:npos;
pos_z=1:npos;

nmid_rel=tmax+2; %relative mid position in the below a arrays which corresponds to zero position in the absolute sense
nstart_rel=0;%relative start position in the below a arrays which corresponds to least -ve position in the absolute sense
nend_rel=0;%relative end position in the below a arrays which corresponds to max +ve position in the absolute sense

a_n_u_prev=zeros(npos,npos,npos); %a's for spin up at each N at previous time
a_n_d_prev=zeros(npos,npos,npos);%a's for spin down at each N at previous time
a_n_u_cur=zeros(npos,npos,npos);%a's for spin up at each N at current time
a_n_d_cur=zeros(npos,npos,npos);%a's for spin down at each N at current time

%%%%%%%%%%%%%%%%%%%%%% SIGMA CODE %%%%%%%%%%%%%%%%%%
%sigma=zeros(1,npos); %the standard deviation of the spatial distribution
sigma=zeros(1,length(time));
sigma_up=zeros(1,length(time));
sigma_down=zeros(1,length(time));

avg_pos=zeros(1,length(time));
avg_pos_up=zeros(1,length(time));
avg_pos_down=zeros(1,length(time));

sig=zeros(1,length(time));
sig_up=zeros(1,length(time));
sig_down=zeros(1,length(time));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% END SIGMA %%%%%%%%%%%%%%%%%%

x_track=[nmid_rel];
y_track=[nmid_rel];
z_track=[nmid_rel];

% T=tmax/2:100:tmax;
T=time;
pos_n_counter=0;
pos_m_counter=0;
pos_l_counter=0;
a=0;
epsilon = 0.1;
p_n_tmax = zeros(npos,npos,npos);
p_n_u_tmax = zeros(npos,npos,npos);
p_n_d_tmax = zeros(npos,npos,npos);

%%%%%%%%%%%%%% INIT POSITIONS %%%%%%%%%%%%%

for posz_iter=1:npos
    for posy_iter=1:npos
        for posx_iter=1:npos
            pos(posz_iter,posy_iter,posx_iter) = sqrt((posx_iter-nmid_rel)^2 + (posy_iter-nmid_rel)^2 + (posz_iter-nmid_rel)^2);
        end;
    end;
end;

%%%%%%%%%%%%%%%%% END INIT POSITIONS %%%%%%%%%%%%%


for iter=1:n_iterations
    display(iter);
    a = 1/(sqrt(2));
%     a = rand;
    a_n_u_prev(nmid_rel,nmid_rel,nmid_rel)=a;
%     a_n_d_prev(nmid_rel,nmid_rel)=(sqrt(1-a^2))*i;
    a_n_d_prev(nmid_rel,nmid_rel,nmid_rel)=a*i;
    
    
    %for all times
    for t=time
        display(t);
        pos_m_counter = 0;
        pos_n_counter = 0;
        nstart_rel = nmid_rel - t;
        nend_rel = nmid_rel + t;
%         if(t ~= 1)
%              a=1/sqrt(2);
             a=(1/(sqrt(2))) + epsilon*rand;
% %             a = epsilon*rand;
%         end;
        randNo = rand;
        if(randNo < 0.333)
            x_track = nstart_rel:nend_rel;
            for pos_l_counter=z_track
                for pos_n_counter=y_track            
                    for pos_m_counter=x_track
    %                     a = rand * 0.5;
                        a_n_u_cur(pos_l_counter,pos_n_counter,pos_m_counter) =  ( sqrt(1-a*a) * a_n_d_prev(pos_l_counter,pos_n_counter,pos_m_counter-1) +  a * a_n_u_prev(pos_l_counter,pos_n_counter,pos_m_counter-1) );
                        a_n_d_cur(pos_l_counter,pos_n_counter,pos_m_counter) =  ( sqrt(1-a*a) * a_n_u_prev(pos_l_counter,pos_n_counter,pos_m_counter+1) -  a * a_n_d_prev(pos_l_counter,pos_n_counter,pos_m_counter+1) );                    
                    end;
                end;
            end;
        elseif(0.333 > randNo && randNo <0.666) 
            y_track = nstart_rel:nend_rel;
            for pos_l_counter=z_track
                for pos_n_counter=y_track            
                    for pos_m_counter=x_track
    %                     a = rand * 0.5;
                        a_n_u_cur(pos_l_counter,pos_n_counter,pos_m_counter) =  ( sqrt(1-a*a) * a_n_d_prev(pos_l_counter,pos_n_counter-1,pos_m_counter) +  a * a_n_u_prev(pos_l_counter,pos_n_counter-1,pos_m_counter) );
                        a_n_d_cur(pos_l_counter,pos_n_counter,pos_m_counter) =  ( sqrt(1-a*a) * a_n_u_prev(pos_l_counter,pos_n_counter+1,pos_m_counter) -  a * a_n_d_prev(pos_l_counter,pos_n_counter+1,pos_m_counter) );                    
                    end;
                end;
            end;
        else
            z_track = nstart_rel:nend_rel;
            for pos_l_counter=z_track
                for pos_n_counter=y_track            
                    for pos_m_counter=x_track
    %                     a = rand * 0.5;
                        a_n_u_cur(pos_l_counter,pos_n_counter,pos_m_counter) =  ( sqrt(1-a*a) * a_n_d_prev(pos_l_counter-1,pos_n_counter,pos_m_counter) +  a * a_n_u_prev(pos_l_counter-1,pos_n_counter,pos_m_counter) );
                        a_n_d_cur(pos_l_counter,pos_n_counter,pos_m_counter) =  ( sqrt(1-a*a) * a_n_u_prev(pos_l_counter+1,pos_n_counter,pos_m_counter) -  a * a_n_d_prev(pos_l_counter+1,pos_n_counter,pos_m_counter) );                    
                    end;
                end;
            end;
        end;
        
        %record the probability densities of up,down and all particles
        p_n_u_tmax = abs(a_n_u_cur).^2;
        p_n_d_tmax = abs(a_n_d_cur).^2; 
        p_n_tmax =  p_n_u_tmax + p_n_d_tmax;   
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%% SIGMA CODE %%%%%%%%%%%%%%%%%%%%%%
        if(iter == n_iterations)
            sigma(t) = ( sum(sum( sum(pos.^2 .* p_n_tmax))) - sum(sum(sum(pos.*p_n_tmax))).^2 );
            sigma_up(t) = (sum(sum( sum((pos).^2 .* p_n_u_tmax))) - sum(sum(sum(pos.*p_n_u_tmax))).^2 );
            sigma_down(t) = (sum(sum( sum((pos).^2 .* p_n_d_tmax))) - sum(sum(sum(pos.*p_n_d_tmax))).^2);


%                 avg_pos(t) = sum((pos-nmid_rel).* p_n_tmax);
%                 avg_pos_up(t) = sum((pos-nmid_rel).* p_n_u_tmax);
%                 avg_pos_down(t) = sum((pos-nmid_rel).* p_n_d_tmax);

        end;
        %%%%%%%%%%%%%%%%%%%%% END SIGMA CODE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%        
     
        
        %set the a arrays of prev time with values of that of the current
        %time and a arrays of current time with zero values
        if(t <tmax)
            a_n_u_prev = a_n_u_cur;
            a_n_d_prev = a_n_d_cur;
            a_n_u_cur = zeros(npos,npos,npos);
            a_n_d_cur = zeros(npos,npos,npos);
        end;
        
        if(iter ~= n_iterations) 
            p_n_tmax = zeros(npos,npos,npos);
            p_n_u_tmax = zeros(npos,npos,npos);
            p_n_d_tmax = zeros(npos,npos,npos);
        end;
        
    end;
   
end;

figure(19)
plot(T,sigma,T,sigma_up,T,sigma_down);
xlabel('T:50');
ylabel('\sigma^2(t) for charge and up-spin,down-spin');
title('\sigma^2(t) for charge and up-spin,down-spin');
grid on;

